PVA of the NBI: Demography - Data Exploration
- Hypothesis: We hypothesize that at the present time the northern bald ibis population can survive without further management and release. We predict that the observed demographic rates will ensure population growth and do not differ between the breeding colonies.
- Study area: Austria, Germany and Italy
- Data: Life History data of the NBI
In this script we created different plots for the description of the NBI population. We created Timelines, Population size per year and Barplots for all individuals and the alive females to show the distribution between the sexes, stages, colonies, raising types and the mortality types.
Data
Save some standard variables
Population size per year
Plot the whole Population size per year
# Calculate the number of living individuals per year
popsize_NBI <- tibble(
"2008" = ifelse(test = NBI$hatch_year <= 2008 & (NBI$death_year >= 2008 | NBI$death_year == 0),
yes = 1,
no = 0) %>%
sum(),
"2009" = ifelse(test = NBI$hatch_year <= 2009 & (NBI$death_year >= 2009 | NBI$death_year == 0),
yes = 1,
no = 0) %>%
sum(),
"2010" = ifelse(test = NBI$hatch_year <= 2010 & (NBI$death_year >= 2010 | NBI$death_year == 0),
yes = 1,
no = 0) %>%
sum(),
"2011" = ifelse(test = NBI$hatch_year <= 2011 & (NBI$death_year >= 2011 | NBI$death_year == 0),
yes = 1,
no = 0) %>%
sum(),
"2012" = ifelse(test = NBI$hatch_year <= 2012 & (NBI$death_year >= 2012 | NBI$death_year == 0),
yes = 1,
no = 0) %>%
sum(),
"2013" = ifelse(test = NBI$hatch_year <= 2013 & (NBI$death_year >= 2013 | NBI$death_year == 0),
yes = 1,
no = 0) %>%
sum(),
"2014" = ifelse(test = NBI$hatch_year <= 2014 & (NBI$death_year >= 2014 | NBI$death_year == 0),
yes = 1,
no = 0) %>%
sum(),
"2015" = ifelse(test = NBI$hatch_year <= 2015 & (NBI$death_year >= 2015 | NBI$death_year == 0),
yes = 1,
no = 0) %>%
sum(),
"2016" = ifelse(test = NBI$hatch_year <= 2016 & (NBI$death_year >= 2016 | NBI$death_year == 0),
yes = 1,
no = 0) %>%
sum(),
"2017" = ifelse(test = NBI$hatch_year <= 2017 & (NBI$death_year >= 2017 | NBI$death_year == 0),
yes = 1,
no = 0) %>%
sum(),
"2018" = ifelse(test = NBI$hatch_year <= 2018 & (NBI$death_year >= 2018 | NBI$death_year == 0),
yes = 1,
no = 0) %>%
sum(),
"2019" = ifelse(test = NBI$hatch_year <= 2019 & (NBI$death_year >= 2019 | NBI$death_year == 0),
yes = 1,
no = 0) %>%
sum())
# Calculate the change per year
for (i in 2:ncol(popsize_NBI)) {
popsize_NBI[2,i] <- popsize_NBI[1,i] - popsize_NBI[1,i-1]
}# only numbers
pop_num <- as.numeric(popsize_NBI[1,])
# plot it
par(mar = (c(5, 4.5, 1, 2) + 0.1))
plot(pop_num, type = "n", xaxt = "n", yaxt = "n", xlab = "Year", ylab = "Number of individuals", cex.lab = 2, ylim = c(0,200))
axis(1, at = 1:12, labels = 2008:2019, cex.axis = 2)
axis(2, cex.axis = 2)
lines(zero_line, col = "grey85")
polygon(c(1:12, rev(1:12)), c(pop_num, rev(zero_line)), col = "grey85", border= NA)
lines(pop_num, lwd = 3)
points(pop_num, pch = 16, cex = 1.5)
box()Save the plot
# png(filename = here("plots", "01_description", "Popsize_per_year.png"))
# par(mar = (c(5, 4.5, 1, 2) + 0.1))
#
# plot(pop_num, type = "n", xaxt = "n", yaxt = "n", xlab = "Year", ylab = "Number of individuals", cex.lab = 2, ylim = c(0,200))
# axis(1, at = 1:12, labels = 2008:2019, cex.axis = 2)
# axis(2, cex.axis = 2)
#
# lines(zero_line, col = "grey85")
# polygon(c(1:12, rev(1:12)), c(pop_num, rev(zero_line)), col = "grey85", border= NA)
#
# lines(pop_num, lwd = 3)
# points(pop_num, pch = 16, cex = 1.5)
#
# box()
# par(mar = mar_default)
#
# dev.off()Customized barplots
Overall population
Show the values for each category and build a matrix with them. The NAs between some values are there for some space between the categories
# Comparison between delivered individuals and all other individuals
NBI_deliv <- NBI[NBI$mortality_cause == "delivered",] # 9 delivered individuals (delivered to a colony outside the project)
NBI_proj <- NBI[NBI$mortality_cause != "delivered",] # 375 not delivered individuals
# Comparison of the reached stages and the life status (alive or dead)
table(NBI_proj$life_status, NBI_proj$reached_stage)##
## 1 2 3 4
## 0 71 28 13 31
## 1 123 45 33 31
##
## 1 2 3 4
## 0 8 1 0 0
## 1 0 0 0 0
##
## 1 2 3 4
## 202 74 46 62
##
## f m u
## 184 195 5
##
## Burghausen Kuchl allocation Ueberlingen Rosegg
## 127 89 57 89 22
##
## foster parent parent supplementation
## 213 162 9
##
## alive delivered during human care electrocution injury/illness other predation shot unknown
## 143 9 4 46 23 2 15 18 124
NBI_mat <- matrix(data = c(
202, 74,46,62,NA, # Stage
184,195, 5,NA, # Sex
127, 89,89,22,57,NA, # Colony (B, K, Ue, R, Allo)
213, 162,9, NA, # Raising type
4,46,23,15,18,2,124,9,143)) # Mortality type (d h c, elec, inj, pred, shot. other, unk, del, alive)Plot
par(mar = c(5,5,2,0))
barplot(NBI_mat, beside = T, ylim = c(0,220), ylab = "Number of Individuals", border = NA,
cex.axis = 3, cex.lab = 3,
col = c(
magma(4, begin = 0.25, end = 0.6), "white", # Stage
plasma(1, direction = -1, begin = 0.0, end = 0.6), viridis(1, begin = 0.2, alpha = 0.9), # Sex
magma(1, begin = 0.1, alpha = 0.6),"white", # Sex
inferno(5,begin =0.65, end = 0.95, direction = -1), "white", # Colony
viridis(3, begin = 0.65, end = 0.85),"white", # Raising type
gray(0:7 / 8), viridis(1, begin = 0.5))) # Mortality type
axis(1, at = c(3,7.5,12.5,17.5,24.5), labels = c("", "", "", "", ""), cex.axis = 3)
axis(1, at = c(3,7.5,12.5,17.5,24.5), labels = c("Stage", "Sex", "Colony", "Raising type", "Mortality type"), cex.axis = 3, line = 1, lwd = 0)
box()Save the Plot
# png(filename = here("plots", "01_description", "Barplot_all.png"), width = 1300, height = 835)
# par(mar = c(5,5,2,0))
#
# barplot(NBI_mat, beside = T, ylim = c(0,220), ylab = "Number of Individuals", border = NA,
# cex.axis = 3, cex.lab = 3,
# col = c(
# magma(4, begin = 0.25, end = 0.6), "white", # Stage
#
# plasma(1, direction = -1, begin = 0.0, end = 0.6), viridis(1, begin = 0.2, alpha = 0.9), # Sex
# magma(1, begin = 0.1, alpha = 0.6),"white", # Sex
#
# inferno(5,begin =0.65, end = 0.95, direction = -1), "white", # Colony
#
# viridis(3, begin = 0.65, end = 0.85),"white", # Raising type
#
# gray(0:7 / 8), viridis(1, begin = 0.5))) # Mortality type
#
# axis(1, at = c(3,7.5,12.5,17.5,24.5), labels = c("", "", "", "", ""), cex.axis = 3)
# axis(1, at = c(3,7.5,12.5,17.5,24.5), labels = c("Stage", "Sex", "Colony", "Raising type", "Mortality type"), cex.axis = 3, line = 1, lwd = 0)
# box()
#
# dev.off()Legend
par(mar = mar_small)
plot(1, type = "n", axes = "F", xlab = "", ylab = "") # empty plot
legend(x = "center",
legend = c("Stage","Juveniles","1 Year Old","2 Year Old","Adults", "",
"Sex","Females","Males","Unknown", "",
"Colony","Burghausen","Kuchl","Überlingen", "Rosegg", "unassigned yet", "",
"Raising type","Foster parents","Biological Parents", "Supplementation", "",
"Mortality type","During human care","Electrocution","Injury or illness","Predation","Shot","Other", "Unknown", "Delivered", "Alive"),
fill = c("white",magma(4, begin = 0.25, end = 0.6), "white", # Stage
"white",plasma(1, direction = -1, begin = 0.0, end = 0.6), viridis(1, begin = 0.2, alpha = 0.9), # Sex
magma(1, begin = 0.1, alpha = 0.6),"white", # Sex
"white",inferno(5,begin =0.65, end = 0.95, direction = -1), "white", # Colony
"white",viridis(3, begin = 0.65, end = 0.85), "white", # Raising type
"white" , gray(0:7 / 8), viridis(1, begin = 0.5)), # Mortality type
border = NA, cex = 3)Save the Legend
# png(filename = here("plots", "01_description", "Barplot_all_legend.png"), width = 2174, height = 1800)
#
# par(mar = mar_small)
# plot(1, type = "n", axes = "F", xlab = "", ylab = "") # empty plot
#
# legend(x = "center",
# legend = c("Stage","Juveniles","1 Year Old","2 Year Old","Adults", "",
# "Sex","Females","Males","Unknown", "",
# "Colony","Burghausen","Kuchl","Überlingen", "Rosegg", "unassigned yet", "",
# "Raising type","Foster parents","Biological Parents", "Supplementation", "",
# "Mortality type","During human care","Electrocution","Injury or illness","Predation","Shot","Other", "Unknown", "Delivered", "Alive"),
# fill = c("white",magma(4, begin = 0.25, end = 0.6), "white", # Stage
# "white",plasma(1, direction = -1, begin = 0.0, end = 0.6), viridis(1, begin = 0.2, alpha = 0.9), # Sex
# magma(1, begin = 0.1, alpha = 0.6),"white", # Sex
# "white",inferno(5,begin =0.65, end = 0.95, direction = -1), "white", # Colony
# "white",viridis(3, begin = 0.65, end = 0.85), "white", # Raising type
# "white" , gray(0:7 / 8), viridis(1, begin = 0.5)), # Mortality type
# border = NA, cex = 3)
#
# par(mar = mar_default)
#
# dev.off()Alive female individuals
##
## 0 1
## f 74 110
## m 69 126
## u 0 5
##
## 1 2 3 4
## 37 11 8 18
##
## Burghausen Kuchl allocation Ueberlingen Rosegg
## 19 15 6 28 6
##
## foster parent parent supplementation
## 41 32 1
NBI_mat_f <- matrix(data = c(
37, 11,8,18,NA, # Stage
19, 15,28,6,6,NA, # Colony (B, K, Ue, R, Allo)
41, 32,1)) # Raising typePlot
par(mar = c(5,5,2,0))
barplot(NBI_mat_f, beside = T, ylim = c(0,50), ylab = "Number of Individuals", border = NA,
cex.axis = 3, cex.lab = 3,
col = c(
magma(4, begin = 0.25, end = 0.6), "white", # Stage
inferno(5,begin =0.65, end = 0.95, direction = -1), "white", # Colony
viridis(3, begin = 0.65, end = 0.85) # Raising type
))
axis(1, at = c(3,8.5,13.5), labels = c("", "", ""), cex.axis = 3)
axis(1, at = c(3,8.5,13.5), labels = c("Stage", "Colony", "Raising type"), cex.axis = 3, line = 1, lwd = 0)
box()Save the Plot
# png(filename = here("plots", "01_description", "Barplot_alive_f.png"), width = 1300, height = 835)
# par(mar = c(5,5,2,0))
#
# barplot(NBI_mat_f, beside = T, ylim = c(0,50), ylab = "Number of Individuals", border = NA,
# cex.axis = 3, cex.lab = 3,
# col = c(
# magma(4, begin = 0.25, end = 0.6), "white", # Stage
#
# inferno(5,begin =0.65, end = 0.95, direction = -1), "white", # Colony
#
# viridis(3, begin = 0.65, end = 0.85) # Raising type
#
# ))
# axis(1, at = c(3,8.5,13.5), labels = c("", "", ""), cex.axis = 3)
# axis(1, at = c(3,8.5,13.5), labels = c("Stage", "Colony", "Raising type"), cex.axis = 3, line = 1, lwd = 0)
#
# box()
#
# dev.off()Legend
par(mar = mar_small)
plot(1, type = "n", axes = "F", xlab = "", ylab = "")
legend("center", legend = c("Stage","Juveniles","1 Year Old","2 Year Old","Adults", "",
"Colony","Burghausen","Kuchl","Überlingen", "Rosegg", "unassigned yet", "",
"Raising type","Foster parents","Biological Parents", "Supplementation"),
fill = c("white",magma(4, begin = 0.25, end = 0.6), "white",
"white",inferno(5,begin =0.65, end = 0.95, direction = -1), "white",
"white",viridis(3, begin = 0.65, end = 0.85)), border = NA, cex = 1)Save the legend
# png(filename = here("plots", "01_description", "Barplot_alive_f_legend.png"), res = 130) # width = 300, height = 400
# par(mar = mar_small)
#
# plot(1, type = "n", axes = "F", xlab = "", ylab = "")
# legend("center", legend = c("Stage","Juveniles","1 Year Old","2 Year Old","Adults", "",
# "Colony","Burghausen","Kuchl","Überlingen", "Rosegg", "unassigned yet", "",
# "Raising type","Foster parents","Biological Parents", "Supplementation"),
# fill = c("white",magma(4, begin = 0.25, end = 0.6), "white",
# "white",inferno(5,begin =0.65, end = 0.95, direction = -1), "white",
# "white",viridis(3, begin = 0.65, end = 0.85)), border = NA, cex = 1)
#
# par(mar = mar_default)
# dev.off()Session Info
## DO NOT REMOVE!
## We store the settings of your computer and the current versions of the
## packages used to allow for reproducibility
Sys.time()## [1] "2022-01-06 19:33:39 CET"
## R version 4.0.3 (2020-10-10)
## Platform: x86_64-w64-mingw32/x64 (64-bit)
## Running under: Windows 10 x64 (build 22000)
##
## Matrix products: default
##
## locale:
## [1] LC_COLLATE=German_Germany.1252 LC_CTYPE=German_Germany.1252 LC_MONETARY=German_Germany.1252 LC_NUMERIC=C LC_TIME=German_Germany.1252
##
## attached base packages:
## [1] splines stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] viridisLite_0.3.0 eeptools_1.2.4 forcats_0.5.0 stringr_1.4.0 dplyr_1.0.2 purrr_0.3.4 readr_1.4.0 tidyr_1.1.2 tibble_3.0.3
## [10] ggplot2_3.3.2 tidyverse_1.3.0 pwr_1.3-0 plot3D_1.4 MuMIn_1.43.17 modelsummary_0.6.5 mgcv_1.8-33 nlme_3.1-149 lattice_0.20-41
## [19] here_0.1 ggeffects_1.0.0 gam_1.20 foreach_1.5.1 effects_4.2-0 carData_3.0-4
##
## loaded via a namespace (and not attached):
## [1] readxl_1.3.1 backports_1.1.10 Hmisc_4.4-1 sp_1.4-4 usethis_1.6.3 digest_0.6.25 htmltools_0.5.0 fansi_0.4.1
## [9] magrittr_1.5 checkmate_2.0.0 memoise_1.1.0 cluster_2.1.0 remotes_2.2.0 modelr_0.1.8 sysfonts_0.8.1 rmdformats_0.3.7
## [17] prettyunits_1.1.1 jpeg_0.1-8.1 colorspace_1.4-1 blob_1.2.1 rvest_0.3.6 mitools_2.4 haven_2.3.1 xfun_0.18
## [25] tcltk_4.0.3 callr_3.5.0 crayon_1.3.4 jsonlite_1.7.1 lme4_1.1-25 zoo_1.8-8 survival_3.2-7 iterators_1.0.13
## [33] glue_1.4.2 kableExtra_1.3.1 gtable_0.3.0 webshot_0.5.2 pkgbuild_1.1.0 abind_1.4-5 scales_1.1.1 DBI_1.1.0
## [41] Rcpp_1.0.5 showtextdb_3.0 htmlTable_2.1.0 foreign_0.8-80 Formula_1.2-4 stats4_4.0.3 survey_4.0 vcd_1.4-8
## [49] htmlwidgets_1.5.2 httr_1.4.2 RColorBrewer_1.1-2 ellipsis_0.3.1 pkgconfig_2.0.3 farver_2.0.3 nnet_7.3-14 dbplyr_1.4.4
## [57] utf8_1.1.4 tidyselect_1.1.0 labeling_0.3 rlang_0.4.7 munsell_0.5.0 cellranger_1.1.0 tools_4.0.3 cli_2.0.2
## [65] generics_0.0.2 sjlabelled_1.1.7 devtools_2.3.2 broom_0.7.4 evaluate_0.14 arm_1.11-2 yaml_2.2.1 tables_0.9.6
## [73] processx_3.4.4 knitr_1.30 fs_1.5.0 showtext_0.9 xml2_1.3.2 compiler_4.0.3 rstudioapi_0.11 png_0.1-7
## [81] testthat_2.3.2 reprex_0.3.0 statmod_1.4.35 stringi_1.5.3 highr_0.8 ps_1.4.0 desc_1.2.0 Matrix_1.2-18
## [89] nloptr_1.2.2.2 vctrs_0.3.4 pillar_1.4.6 lifecycle_0.2.0 lmtest_0.9-38 estimability_1.3 maptools_1.0-2 data.table_1.13.2
## [97] insight_0.11.1 R6_2.4.1 latticeExtra_0.6-29 bookdown_0.20 gridExtra_2.3 sessioninfo_1.1.1 codetools_0.2-16 boot_1.3-25
## [105] MASS_7.3-53 assertthat_0.2.1 pkgload_1.1.0 rprojroot_1.3-2 withr_2.3.0 d6_0.1.0.0 hms_0.5.3 grid_4.0.3
## [113] rpart_4.1-15 coda_0.19-4 minqa_1.2.4 rmarkdown_2.4 misc3d_0.9-0 git2r_0.27.1 lubridate_1.7.9 base64enc_0.1-3
## [121] tinytex_0.26